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Factorizable upwind schemes: 
the triangular unstructured grid formulation 

David Sidilkover* 1 

ICASE, NASA Langley Research Center, Hampton, VA 23681 
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The upwind factorizable schemes for the equations of fluid was introduced recently. 

They facilitate achieving the Textbook Multigrid Efficiency (TME) and are expected also 
to result in the solvers of unparalleled robustness. The approach itself is very general. 
Therefore, it may well become a general framework for the large-scale Computational 
Fluid Dynamics. In this paper we outline the triangular grid formulation of the factor- 
izable schemes. The derivation is based on the fact that the factorizable schemes can 
be expressed entirely using vector notation, without explicitly mentioning a particular 
coordinate frame. We describe the resulting discrete scheme in detail and present some 
computational results verifying the basic properties of the scheme/solver. 


Introduction 

This work is a part of effort, going on at NASA 
Langley for several years towards constructing a new 
generat ion of the flow solvers (see, for instance Thomas 
et alb, Roberts et al 2 ). The key idea, that was sug- 
gested by Brandt, 3 is to use a special relaxation that 
recognizes the mixed character of a system of PDEs. 
Then each sub-factor of the system can be treated in 
different (optimal for it) way. Tt is well-known that 
the Euler equations “consist" of two different factors: 
advection and full- potential operators. The advection 
part can be treated very efficiently, sav, by the march- 
ing relaxation. The full-potential operator is of the 
elliptic type in the subsonic regime. Therefore, it can 
be treated very efficiently by multigrid. In the su- 
personic regime it becomes hyperbolic: wave equation 
with respect to the flow direction. For Mach number 
substantially larger than one, the entire system can 
be solved efficiently by marching. Nearly sonic speed 
regime can be dealt with by multigrid, but requires 
some special care. 

Solvers based on such a special (Distributive) relax- 
ation were constructed initially for incompressible flow 
and were based on the staggered grid discretizations. 
The optimal multigrid efficiency was demonstrated 
Brandt and Yavneh. 4 Staggered-grid discretizations 
exist also for the compressible flow equations. How- 
ever, they all are limited to the subsonic regime, since 
they have no shock-capturing capabilities. The stan- 
dard shock-capturing schemes, on the other hand, are 
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not factorizable, i.c. they do not reflect the mixed 
character of the PDEs (see Sidilkover 5 ). Therefore, it 
is not possible to construct, a relaxation of Dist ributive 
type that can be used with these schemes. 

Clearly, there is a need for discretizations that are 
both factorizable and have shock-capturing capabili- 
ties. A factorizable upwind scheme w r as constructed in 
Sidilkover 6 for the case of Cartesian grids. A detailed 
description of its extension to the case of structured 
body-fitted grids is given in Sidilkover et al.' A set of 
numerical results was presented in Roberts et al. 8 

The purpose of this paper is to present a construc- 
tion of a factorizable scheme on triangular unstruc- 
tured grids. 

Euler equations and their properties 

The non-conservative quasilinear formulation of the 
compressible Euler equations in three dimensions can 
be written using vector notat ion as follows 

u ■ Ys - 0 (la) 

pit ■ Vw + Vp = 0 ( lb) 

pc 2 V u + t7 * Vp = 0, (lc) 

where s denotes the entropy and 

ds = dp- -L (2) 

C“ 

it is the velocity vector, the pressure p is given by 

P= (7- ( 3 ) 

the speed of sound 
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In this section wc recall some of the basic properties 
of the equations (1). It is sufficient for the purpose 
of the analysis to assume constancy of the coefficients. 
It is known that this system of equations is of the 
mixed type: it consists of the advection and the Full- 
Potential factors. This can be made obvious by intro- 
ducing the new set of variables. 

Recall that a vector field can be decomposed into 
solenoidal and irrotational parts 

u - V x 0 T (5) 

where <p is the potential and V’ is the streamfunction. 
The pressure gradient is related to the gradient of the 
potential as follows 

dp ~ — pi 1 * Vdf> (6) 

Substituting the new variables (j> and V into the pres- 
sure equation we obtain the Full-Potential equation 

p[r 2 V 2 -(ff- V) 2 ]0 = 0 (7) 

Note, that all the terms involving the streamfunction 

cancel out. 

Performing the variable substitution in the momen- 
tum equations gives 

puV{Vx if) = 0 (8) 

Note, that all the terms involving the potential vari- 

able cancel out. 

Introducing a new variable vorticity 

Q = Af (9) 

and applying operator Vx to (8), we obtain 

puV(l = 0 (10) 

This verifies indeed that the Euler system is of the 
mixed type. The advection factor is represented by 
the equations for entropy (la) and vorticity (10). The 
full-potential factor is given by (7). It also makes it 
clear that (for the linear constant coefficients case) the 
momentum equations (lb) drive the solenoidal part of 
the solution, while the irrotational part of the solution 
is subject solely to the pressure equation (lc). 

In a general nonlinear case (away from singularities, 
like shocks and contact discontinuities) there is a weak 
coupling between different factors due to the so-called 
subprincipal terms. This coupling can be neglected 
for the purpose of the construction of a fast solver (see 
Brandt 3 ). Therefore, the latter can rely entirely on 
the analysis of the linear case. 

Preparations for the scheme 
construction 

When constructing a discrete approximation to the 
Euler equations, the central scheme can serve a, ba- 
sic building block. However, it is crucial to include 


a certain artificial dissipation in the discretization for 
stability reasons. One of the additional problems than 
becomes how to compensate for the loss of accuracy 
due to the artificial dissipation. Various ways to deal 
with this issue received an extensive coverage in the 
literature. Constructing a factorizable scheme implies 
resolving this issue in a very specific way. 

FDA analysis 

The First Differentia] Approximation (FDA) (or the 
modified equations) corresponding to a certain discrete 
scheme is the PDFs augmented by the leading error 
terms. 

We shall start our analysis with formulating the 
FDA for the factorizable genuinely multidimensional 
scheme. The observation made in regarding the gen- 
uinely multidimensional upwind scheme introduced in 
Sidilkover 9 was that a part of the artificial dissipa- 
tion present in the discrete momentum equations (in 
subsonic case) is proportional to the gradient of the 
residual of the pressure equation. The artificial dissi- 
pation of the pressure equation is proportional to the 
divergence of the residuals of the momentum equa- 
tions. A vector formulation of the entire scheme (on 
the Cartesian grids) is given in Sidilkover. 10 

The fact that the entire scheme can be expressed 
using the vector notation appeared to be very instru- 
mental for the purpose of extending the factorizable 
to the structured body-fitted grids (See Sidilkover et 
al 7 and Roberts et al 8 ). It is of very important for the 
purpose of this paper too. 

The FDA of a factorizable scheme for the Euler 
equations is given by the following 

q.s = 0 (11 a) 

pqii +Vp- — - V ■ u + u • Vp) - VD 
2 c 

= 0 (lib) 

pe 2 V • u + u ■ Vp — ^|^cV * (pd ■ Vu -f Vp) 

= 0 (lie) 

The term VD in the momentum equations plays an 
important role when the operator q is discretized us- 
ing an advection scheme of a certain type to maintain 
the second order accuracy and factorizability of the 
whole Euler scheme. This special type of the advection 
scheme allows to upgrade the accuracy of the advec- 
tion factor to the second order without affecting the 
discrete full- potential part. For now we omit the term 
V D and consider q to correspond to the standard first 
order accurate advection scheme. 

A factorizable scheme corresponding to the FDA as 
given by (11) is stable for subsonic case only. The 
scheme can be extended so it will be valid for tran- 
sonic/supersoni c regime by a simple modification: in- 
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troduction of a cut-off parameter on the terms involv- 
ing pressure in the artificial dissipation 

q s — 0(12a) 

pqu -f Vp - -~-V{pc~V * t7+ - u ■ Vp) -VD 

2 c K 

= 0(12b) 

pc 2 V * i! + u ■ Vp ^j-cV * Vu 4- K*Vp) 

= 0, (12c) 

where 

k = max( 1 , A/ 2 ) (13) 

We shall demonstrate the factorizability property of 
FDA corresponding to the new scheme. Introduce the 
auxiliary variables - potential <j> and streamfunction i- 
and substitute the following expressions for the vari- 
ables *7, p 

u = (I- ^--u- V)V X if 
2 c 

(14a) 

2 c 

dp = -p(u ■ V - ^-cV 2 )# (14b) 

into (II). Introducing the vorticity variable 

n = Vx{! _ x f (15) 

and applying the V x operator to the momentum equa- 
tions we obtain 

p?7 ■ VQ = 0. (16) 

Note, that, as well as in the PDE case, all the terms 
involving the potential variable canceled out. Substi- 
tuting the auxiliary variables into the pressure equa- 
tion, it is easy to verify that all the terms involving 
the streamfunction cancel. 

We can summarize that due to the specific form of 
the artificial dissipation, the FDA (11) of the discrete 
scheme is factorizable - it reflects the mixed character 
of the original system of PDEs. 

Special discrete operators 

We consider the two-dimensional case from now on. 
When discretizing the derivatives, a special care needs 
to be taken of what kind of discrete operators are used 
in order to preserve the factorizability property at the 
discrete level. 

Note, that when demonstrating the factorizability 
property of the scheme’s FDA (11) we used the facts 
of the following type 

d xx d y - d xy d x 

dyydjr = drydy ( 17 ) 

0 X y0 X y 0 XX 0yy 


In order to obtain the factorizable discrete scheme, 
we need to introduce some finite differences that pos- 
sess the property analogous to (17). Such finite differ- 
ences were introduced for the case of structured grids 
in. 6 

We have to find such differences for the case of 
triangular grids (see Figl), where (a\ t/) is a local 
(noil-orthogonal) frame. We shall illustrate this on 
a simple example. Consider approximating the par- 
tial derivative d xx . Assuming that we have at our 
disposal the '•compact'" stencil that involves 7 points: 
0, 1,2, 3, 4, 5, 6, there is only one way of doing it, 
namely by d xr defined as follows 

Q xx u - ( w 4 - 2«o + mi )/fr ( 1 8) 

Differences defined in such a way do not have the 
property (17). We can conclude that using the com- 
pact stencil 7-point stencil only there is no way to 
achieve this property. We know from 6 that the 9-point, 
box structured grid stencil is sufficient for this pur- 
pose. There are several ways to augment the compact 
7-point stencil to the 9-point one in the current trian- 
gular grid context . It is clear that the 7-point stencil, 
therefore, needs to be augmented. We can do it by- 
adding to it. 6 more nodes: 7,8,9, 10, 11, 12. The par- 
tial derivative d XT can then be approximated by a wide 
difference 

B xx — [(«4 — 2 My + m . i )/2 

+ (u-io + « 3 - u 2 - us)/ 8 (19) 

4 (u u 4 «6 “ - w 7 )/8]//r 

The derivatives d yy , d x . d y can be approximated in the 
analogous way. It is easy to verify that such differences 
possess the property (17). 

Structure of some artificial dissipation terms 

The central part of the scheme is the constructed 
in the standard way. The artificial dissipation terms 
corresponding to the advection scheme are evaluated 
in the standard fashion as well. A special care needs 
to be taken of the other artificial dissipation terms. 

Recall, that the artificial dissipation terms in the 
momentum equations (lib) that are subject to the 
gradient operator are the residual of the pressure equa- 
tion. Denote them 

R p — pc 2 V ■ u + ii • Vp (20) 

and the expression subject to the action of the diver- 
gence operator in the pressure equation (lie) is the 
residual of the momentum equations 

R m = pit • Vu + Vp (21) 
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Fig, 1 Computational grid segment and a control volume. 


Constructing the discrete scheme 

Our goal now is to derive a discretization of a con- 
servation law (scalar and a system). For this purpose 
we need to evaluate the numerical fluxes through each 
of the faces of the dual-median cell (see Fig.l). 

Global and local coordinate frames 



where (o^,/^) and (a rn 3 }] ) are the unit vectors in 
the direction of the £ and // coordinate axes respec- 
tively. The relationship between the Cartesian and 
contravariant velocity components is described as fol- 
lows 



The Jacobian of this coordinate rotation 

J = det H = a^ n - faa n (25) 



Fig. 2 Nodes used for evaluating the flux through 
a face of dual-median cell. 
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The inverse of the Hessian H 


H 



A» -®x ^ 
q c y 


(26) 


It is convenient to use the scaled contravariant velocity 
components 


U — Ju — u3 1} - vq 7) 
V = Jv = - uj3^ 4 


Th relationship between the Cartesian and covariant 
velocity components is given by the following 


H t 


u 


v 



(28) 


U — ua £ + vs3^ 
V = ua n 4- vflq 


(29) 


The covariant and covariant velocities are related as 
follows 

« T " (*)=(") m 

H t H=( , + \ (31) 

\ a ^°r/ 4 P$P T} Ck r] 4 Pf) ) 

The total velocity squared 

\u\ 2 — u 2 4 v 2 = uU 4 vV (32) 

Scalar advection 

Consider a scalar advection equation 

s t + ns x 4 vs y = 0. (33) 

The discrete equation to solve for s at point 0 is ob- 
tained by balancing fluxes through the surface of the 
dual median cell (see Fig-1) 

iVj 

£[/ h„]i = 0, (34) 

» = 1 

where h v is the length of the corresponding face and 
the numerical flux 

f — — 4 2 ^* ( s ° + tS ?) ( 35 ) 

where d £ stands for a divided difference 

h^d^s = s 2 - s 0 (36) 

The discrete equation to solve for s at point 0 is ob- 
tained by substituting numerical fluxes evaluated by 
analogy to (35) on all the faces of the dual- median cell 
and substituting them into the flux-balance equation 
(34). 


Euler system 

Integrating the Euler equations in the conservative 
form over the control- volume (dual-median cell) and 
applying the Green's theorem, we obtain 



F • ndl = 0, 


(37) 


where dC is the control volume’s boundary n is a unit 
vector normal to the boundary and 


. 



( pv \ 

F = 

pir + p 

i 4 

puv 
*> , 


puv 


pv- + p 


(E +p)u ) 


(E + p)v ) 


(38) 


A conservative discretization of the Euler system can 
be written in the following form 


= °- ( 39 ) 

i=i 

The numerical flux through a face can be represented 
as a sum of central and the artificial dissipation parts 

F= F c + F d (40) 


The central portion of the flux is given by 

F* c = i[ F(uj.) + F(u*)]-»7, (41) 

where now n — — a, ? ) is a unit vector normal to 

the face. The rest of this section is dedicated to the 
question of deriving the diffusive portion of the numer- 
ical fluxes. 

We would like to emphasize that it is fairly simple 
to implement the new discretization within the exist- 
ing control- volume computer codes. It requires only 
the new numerical flux routine. Such a routine can be 
written in several simple steps, starting from the stan- 
dard upwind scheme and performing the modifications 
gradually. 


The standard upwind scheme 
The first step is to rewrite the standard upwind for 
the subsonic case without explicit ment ion of the char- 
acteristic variables, etc. The scheme is given by the 
following numerical flux 



h \ 

f i 

h 

/i / 


/ \C\P* \ 

1 pcd*U + U/cd* P 

~ 2 H _r\c\d£v 

\ pcUdfU + cdjlp / 


(42) 


The derivation presumes that these fluxes correspond 
to the local orthogonal coordinate frame associated 
with the cell-face. Therefore, the momentum equa- 
tions diffusive fluxes can be rotated to the global co- 
ordinate frame (.r, y) by in the following wav 
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The obtain diffusive fluxes correspond to the quasi- 
linear nonconservative formulation of the Euler equa- 
tions, The need to be transformed, therefore, to the 
form appropriate for the conservative discretization 

F d = Mf d (44) 

where 


1 

0 

0 

1/e 2 

\ 

u 

1 

0 

u/c 2 


V 

0 

1 

v/c 2 


(« ■ «)/ 2 

u 

V 

1 /(7- 1) + (« • u)/(2c 2 ) 

/ 


(45) 


A modified scheme 

As an intermediate step towards constructing the 
factorizable scheme we can consider the following case: 



h \ 


\p\ Qs \ 

f d = 

h 

h 

=-;*« 

pcd'pr + U/cd^p 
p\V\%V 


/•i ! 


pcUqU + e%p J 


The main difference from the standard scheme is that 
the momentum equations diffusive fluxes (the second 
and third components) are now attributed to the mo- 
mentum equations in the covariant directions (£ and 
77 ). Therefore, the transformation back to the global 
orthogonal frame takes the following form 

( h ) = , ' ,Tr ' ( 1 ) <47 > 

This change is necessary towards eventually obtaining 
the approximation for the gradient operator term in 
(lib) 

Another intermediate step 



h \ 

h 

h 

h / 


/ i M* \ 

1 pedj! U + U / cdj: p 

~2* p\U\d*V 

\ pcUd'u + jcd^p / 


(48) 


The coefficient in front of the pressure difference in 
the fourth component constitutes the change from the 
previous case. It is necessary in order arrive later to 
the approximation of the divergence term 


The factorizable scheme 

Now we shall incorporate the correction (mixed 
derivative) terms into the scheme, so that some of 
the terms in the artificial dissipation can be viewed as 
residuals of the momentum and the pressure equations 
and, therefore, are second order small. The latter mod- 
ification together with introducing the wide differences 
that possess the commutativity properties results in a 


factorizable scheme. Implementation of all these steps 
can be done in several steps as well, testing the routine 
after each modification. 

We start from re- interpreting of some standard no- 
tions. A st andard “narrow 51 divided difference can be 
defined also as 


^ S 012 ^ 0l2 + S° 23 3f 3 

dc — 


£° 12 +5 i 


'023 


(49) 


Define a “wide divided difference” 

d } ) = (2 S 012 ch 012 + 2 S 023 3f 3 

" +S° 61 dJ 61 + S I82 ^ & 

+5 293 c)f 3 + 5 034 4 34 )/ 

(2S 012 + 2S 023 + £P 61 -f S 182 + S 293 + S 034 ) 

(50) 


Similarly, we can define Off and . Adding the some 
specific ^/-derivative terms (all the differences used here 
are narrow) to the artificial dissipation of the pressure 
equation obtain the residual of the momentum equa- 
tion in the direction normal to the cell face 

Rh = />(U 3 f + V^)U + + ( Q ( a » + 

(51) 

We also add some ^-derivative terms to the artificial 
dissipation of the momentum equat ion in £ direction 
to obtain the residual of the pressure equation. All 
the differences used here are wide. The notation (IR 
for the residual) reflects this fact 

IR* = pc 2 (%U + B h v V) + (Ud'l + Vd h n )p (52) 
The artificial dissipation then takes the following form 


h \ 



h 

1 

+ h iP \V\%a) 

h 

- ~2 

h(P\U\$V 

\ut 

<T P IC Rfc / 


(53) 


The need for rescaling the artificial dissipation in order 
to avoid the quasi-ellipticity of the full-potential fac- 
tors approximation for the low-speed flow regions was 
established in . 6 The scaling parameters <r rn , a p serve 
this purpose. In the conducted preliminary numerical 
test (see next section) they were taken to be 1. The 
parameter / was taken to be equal to h 

There remains a need to introduce some modifica- 
tions into the central part of the scheme. It is neces- 
sary for factorizabihty that the pressure gradient term 
in the momentum equations is approximated by dif- 
ferences. Using wide differences to approximate the 
pressure advection operator in the pressure equation 
is not necessary for factorizabihty but is still benefi- 
cial since it results in a better form of the discrete 
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Fig. 3 The computational grid. 



Fig. 4 Solution: density contours. 


full-potential factor. These modifications can be in- 
troduced in a very simple way using a trick by Tom 
Roberts: adding certain terms to the artificial dissipa- 
tion. Introduce the following undivided difference 

s 1 ' = [(2S 023 cL 023 + s 293 d 293 + s° 34 d° v 3i ) 

—(2 S° 12 0“ is + S 061 d ° 61 + 5 182 ^ 82 )]/ 

2 hf 


(54) 


The artificial dissipation terms are then augmented 
a.s follows 


h 

Sa 


h + §™P 

/-j + 


(55) 


Returning to the global Cartesian coordinate frame 

(56) 


Wr 1 


h 

h 


and converting to the conservative form 
F d = Mf d 


(57) 


This describes the scheme that was used in the pre- 
liminary numerical experiment reported below. 


Preliminary numerical results 

The work on implementing the new scheme within 
the FUN2D 11 code has just began. Our very first aim 
is just to implement the numerical fluxes and to verify 
the correctness of the residual evaluation. 

The testcase presented is a subsonic flow (Mach = 
.2) in a channel with a bump. The grid consists of 
1375 nodes (see Fig. 3). A second order version of the 
new scheme was used. The contour plots of density 
are presented in Fig. 4. We also present for compar- 
ison in Fig. 5 density contours of the solution to the 
same problem using the standard second order scheme 
with limiters switched off. order upwind scheme. The 
solution obtained using the new scheme is at least as 
accurate as the one using the standard scheme. 


Fig. 5 Solution obtained using the second order 
upwind scheme: density contours. 
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